Stress–Charge Nonlinear Physical Description and Tensor Symmetries for Piezoelectric Materials

Nonlinear piezoelectric materials are raised as a great replacement for devices that require low power consumption, high sensitivity, and accurate transduction, fitting with the demanding requirements of new technologies such as the Fifth-Generation of telecommunications (5G), the Internet of Things (IoT), and modern radio frequency (RF) applications. In this work, the state equations that correctly predict the nonlinear piezoelectric phenomena observed experimentally are presented. Furthermore, we developed a fast methodology to implement the state equations in the main FEM simulation software, allowing an easy design and characterization of this type of device, as the symmetry structures for high-order tensors are shown and explained. The operation regime of each high-order tensor is discussed and connected with the main nonlinear phenomena reported in the literature. Finally, to demonstrate our theoretical deductions, we used the experimental measurements, which presented the nonlinear effects, which were reproduced through simulations, obtaining maximum percent errors for the effective elasticity constants, relative effective permittivity, and resonance frequencies of 0.79%, 2.9%, and 0.3%, respectively, giving a proof of the potential of the nonlinear state equations presented for the unifying of all nonlinear phenomena observed in the piezoelectric devices.


Introduction
Piezoelectric materials have been used in several application fields because their performance and set of physical properties meet the requirements in a wide scope of applications. Since the discovery of the piezoelectric effect by the Curie brothers in the 1880s, these types of materials were mainly used in transduction applications, until the 1970s, when their implementation in radio frequency (RF) applications was developed [1], and currently, the semiconductor manufacturing process allows their use in applications where the transduction between mechanical and electric fields is mandatory at the micro-scale.
Amorphous piezoelectric materials are used in applications where miniaturization is not required, and for this reason, currently, crystalline piezoelectric materials dominate the market and industry, mainly with microelectromechanical system (MEMS) devices due to the reproducibility of their physical and system properties. Consequently, the research on these materials is focused on crystalline composites that have high chemical resistance, high breakdown voltages, and high rigidity for RF applications.
Since its discovery, several fabrication techniques have been developed to obtain piezoelectric materials, where the chemical-based techniques have been of interest due to the requirements of thin-film technologies [2]. Deposition techniques such as metal-oxide chemical vapor deposition (MOCVD) [3] and chemical solution deposition (CSD) [4] are current research topics. Furthermore, there are CMOS-compatible deposition techniques, since these processes have low fabrication temperatures, such as sputtering-based techniques, which can obtain high levels of crystallinity [5,6], being an ideal fabrication process to apply the nonlinear phenomena of piezoelectric materials in a new scope of applications [7][8][9][10].
Currently, the main applications of piezoelectric materials are embedded in the MEMS scope, because they use the accurate transduction capability to implement them in several types of applications such as micro and nano-resonators [11,12], energy harvesters [13], accelerometers [14], wearable devices [15], micro-and nano-actuators [16], and sensors for gasses [17] and electrostatic charge [18]. In general, the applications cited share demanding requirements such as low power consumption, high sensitivity, accurate transduction, great chemical resistance, and good enough electrical and mechanical properties, where all of these conditions are met by piezoelectric materials. The modeling of the mentioned devices using the linear description of traditional state equations [19] gives acceptable errors by its predictions; nevertheless, under relatively high electric fields (>10 6 V/m) and deformation, the physical behavior of materials is not predicted correctly [7,20], and the need for a complete first-principles physical description of the nonlinear phenomena for piezoelectric materials emerges as a mandatory tool for new designs in demanding applications of the industry, such as the Fifth-Generation of telecommunications (5G) and the Internet of Things (IoT).
There are applications that use the nonlinear properties with the same targets as the linear applications exposed above such as actuators [21,22], energy harvesters [23], sensors [24], memories [25], and tunable devices [7,26,27]. In all of these works, the physical and electrical behavior of the system is explained through mathematical models [10,[28][29][30][31] or first-principles deductions (a specific thermodynamic formulation) [32][33][34][35], where the models are only valid for a specific geometry disposition or layer stack, while the physical formulations are general, but very difficult to solve analytically. The case of the hysteresis nonlinear effect is a special topic since its behavior has remnant fields after time; its formulation in deformation-charge form and micro-mechanical modeling was exposed in [36,37] respectively. In the models cited, the core concept used is the algebraic or complex expansion of the material parameters, resulting in adjustments of the macroscopic magnitudes of the physicalsystem, e.g., resonance frequency, effective material constants, lumped elements of equivalent circuits, and quality factor, among others. All of these results produce an imbalance of the state equations, being the core problem of models for nonlinear piezoelectric applications since the introduction of adjustment parameters in the material constants reproduces the macroscopic behavior of the phenomena; nevertheless, the physical behavior of the effect is not described by the equations. In contrast, the first-principles formulations are based on the balance of the microscopic states of the physicalsystem, resulting in the prediction of the physical behavior of the macroscopic states, being a complete physical description of the nonlinear effects where the state equations remain balanced; consequently, the solver's calculus is more difficult and time consuming. Due to this, to include this type of device within the integrated circuit (IC) industry, mathematical tools are needed that allow fast design and manufacturing processes, such as the finite-element-method (FEM)-based design accompanied by compatibility with the main IC fabrication processes, such as CMOS, PD-SOI, FD-SOI, and FinFET. In summary, a thermodynamic formulation with easy implementation in leading FEM simulation software (e.g., COMSOL Multiphysics and COVENTOR) is mandatory for the inclusion of nonlinear piezoelectric devices within the semiconductor industry, and that set of demanding characteristics for the mathematical and physics tools are contained in the formulation presented in this work.
To simulate the nonlinearities and physical behavior of the piezoelectric materials, it is necessary to know the nonlinear state equations with an easy methodology to include them in the FEM solver' calculus; consequently, the symmetry structure of high-order tensors must be given as well. Despite this, the methodologies found in the literature to implement nonlinear behaviors in leading FEM simulation software are complicated to carry out, and at the same time, the symmetry structures cannot be found (only some components for a few types of materials [32]). For these reasons, the nonlinear applications reviewed cannot be explained by a unified set of equations with known tensor structures, making the industrial adoption of these types of applications more difficult despite their advantages.
Taking into account the above discussion, in this work, we present a complete physical description of the nonlinear behavior of piezoelectric materials, obtained through the deduction from first-principles of the nonlinear state Equations (until third-order phenomena), the transformation laws required, and the symmetry structures of the tensors, for each of the thirty-two point groups of symmetry (all types of crystalline materials). Furthermore, a methodology with an easy way to implement the state equations and highorder tensors components in the main FEM simulation software is presented, allowing designing and manufacturing devices that can be used in the 5G, IoT, and RF application scopes. Finally, this work gives the MEMS scientific community all the mathematical and physics tools needed to research new types of applications and optimizations for nonlinear piezoelectric devices.

Stress-Charge Nonlinear Formulation
A suitable thermodynamic representation for including the nonlinear effect within FEM simulators is the stress-charge formulation due to the characteristics of direct solvers, since the physical behavior of the electrical permittivity and elasticity constants are wellknown parameters of crystalline materials; in the literature can be also found references to perform the energy and dissipation calculus [38].
The following deductions are focused on crystalline materials. The theoretical development starts from first-principles using the Voigt form for mechanical tensors, the Einstein sum convention,and the recommended notation for point groups of symmetry by the International Union of Crystallography (IUCr) [39]. The entropy and the temperature contributions were neglected due to the solid phase of materials, the low power dissipation (around 10 mW/mm 2 ), and the nonlinear perturbative operation regime of the devices. In the next sections, we discuss the experimental limits that govern the theoretical development presented.
From the eight possible formulations [19], we used the thermodynamic potential of the electric Gibbs function [40], the total differential of which is defined for the piezoelectric effect as where D k , E k , T λ , and S λ are the electric displacement vector, electric field, stress field, and deformation field, respectively. Therefore, considering the properties of the total differential of a multivariable function, the total differentials for dependent variables are where e T kλ are the piezoelectric coefficients, C λµ the elastic constants, and ij the electrical permittivity.
To deduce the nonlinear formulation, we expanded the tensor coefficients in Equation (2) through a Taylor series centered at zero and took into account that the dependent variables are a function of S λ and E k , then the elastic constants are where the high-order derivatives were measured at constant deformation and the electric field equals zero. Through an analogous procedure, we can obtain all coefficient tensors of Equation (2) as a function of S λ and E k . Since Equation (1) is a total differential, we have then, considering that G 2 is a physical magnitude, it is continuous, has an exact differential, and has derivatives up to third-order, and knowing the mixed derivatives equivalence, we obtain Applying Equations (4) and (5), we can define g λµk as Considering the other tensors' coefficients in Equation (2) and applying the same procedure for Equations (3) to (6), we define the remaining high-order tensors as and Subsequently, by replacing Equations (3), (6), (7), and (8) in (2), we obtain After integrating Equation (9), we finally obtain the nonlinear state equations for the piezoelectric effect considering effects up to third-order In Equation (10), given the equivalence between the Voigt and traditional mechanical notation, an algebraic factor is not needed; this means Equation (10) describes how the exchange of the coupling fields' magnitudes is performed through the direct and converse piezoelectric effect, while the nonlinear contributions generated by the relatively high electric and deformation fields were considered. These conditions subject the material to mechanical and electrical stress, producing variations in all material parameters, as will be shown in the nonlinear effects section.
In this context, t λµν is the contribution to the stress field due to strong deformations taking importance in the plastic operation regime. Furthermore, t λµν relates the orthogonal deformations S µ and S ν that produce the change of the stress field with respect to the linear approximation. Analogously, r ijk is a correction term for the electric displacement vector as function of very high electric fields E j and E k , so this tensor governs the dielectric polarization when |E i | is around 10 9 V/m. g λµk is responsible for the elasto-electric effect (in the literature, also known as nonlinear electrostriction and the electro-elastic effect), since its contribution to the stress field g λµk S µ E k provokes an augmentation of the effective elasticity constants, producing a stiffening of the material. In the same way, q ijλ E j S λ contributes to the electric permittivity due to the strains S λ , and it is responsible for the change in the effective permittivity of a material subject to relatively high electric fields. Finally, The last quadratic terms of Equation (10) are a contribution to the stress and electric displacement field, modifying the value of the coupling piezoelectric coefficients e kλ and e T iµ , respectively. With this approach, the state equations presented remain balanced, while considering the nonlinear effects, and therefore, the physical behavior of the microscopic and macroscopic states of the physics system are predicted correctly.

Transformation Laws
To obtain the symmetry structure of any tensor, we need to know the transformation laws and the symmetry generators a ij of each crystal type (point group of symmetry). Taking into account the recommended notation for crystal classes and point groups by IUCr [39], a ij belongs to special orthonormal group SO(3), since it represents a generic 3D rotation. Furthermore, the transformation laws for the high-order tensors must meet the constraints of the positive energy and generate stable states for the system (e.g., the vanishing of the total torque about the origin), and their symmetry structure must only depend on the point group of the material. Then, to deduce the transformation laws for each high-order tensor, we start with the example of the calculus of the transformation law for the electrical permittivity of any material. The transformations laws for the electric displacement vector and electric field are where the superscript means a transformed magnitude. Now, knowing the law for the polarization of a material: the target is to obtain an equivalent equation in terms of transformed magnitudes, so using Equation (13), we obtain where the transformation law for electrical permittivity is deduced from the symmetry condition, which means that, after transformation, the tensor form (structure) remains invariant: For the q ijλ high-order tensor, we need the transformation law for the deformation field: where N µν is a function of the symmetry generator a ij [38]. Then, using we obtain where the transformation law obtained for the q ijλ tensor is Through an analogous deduction, the transformation laws for the nonlinear tensors in the Equation (10) can be obtained, and they are shown below: At this point, the reader can notice that there are two ways to obtain the transformation laws, one per state equation in (10). Both ways have equivalent results knowing the properties of a symmetry generator a ij (belongs to the SO(3) group) and the M and N matrices [38]:

Symmetry Structure for High-Order Tensors
The structure of the tensors can be calculated from Equation (21), the generator symmetry a ij for the specific crystal type, and a last mathematical constraint based on the mixed derivatives theorem. Then, selecting a transformation law from Equation (21) for the desired high-order tensor structure, a specific point group of symmetry (e.g., 6 mm), and applying Equation (23) in the transformation law selected, we obtain an undetermined algebraic linear system, which, after being solved, we obtain the structure of the tensor in terms of a few unique components, which represents the contribution of the specific tensor to the nonlinear behavior of the piezoelectric material. This procedure to obtain the symmetry structure of the high-order tensors was tested through obtaining the symmetry structure of known tensors for the thirty-two point groups; specifically, the elasticity constants, electrical permittivity, piezoelectric coupling coefficients, and r ijk tensor were reproduced; the last one is the only high-order tensor, whose complete symmetry structure has been published [41]. Table 1 presents a first approximation of the high-order tensors for two common piezoelectrics, PZT − 5H and aluminum nitride (AlN), which belong to the 4 mm and 6 mm point groups, respectively. These results were obtained after reviewing the literature and noticing that, when an excitation signal provokes the appearance of the nonlinear effects [7,10,33], we suppose a variation around 2% for dependent variables with respect to the linear approximation. Table 2 presents the symbols and particular numeration for the thirty-two point groups of symmetry; this numeration is used in the tables where the tensor structures are shown, and all high-order tensor components are introduced only by subscripts. The symmetry structures of tensors q jkλ and g λνi are shown in Tables 3 and 4, respectively; the first column contains the component of the high-order tensor and the following columns its corresponding value for a specific point group. The symmetry structure of q jkλ depends only on the Laue symmetry group. All types of crystalline materials have q jkλ and t λµν different from zero in at least one component, and the g λνi and r ijk tensors are null if the material does not exhibit linear piezoelectric behavior (this means they are centrosymmetric crystals), with the only exception of point group 432, where g λνi is not zero and r ijk remains null. The symmetry structure of t λµν for some point groups is shown in Appendix A in Table A1; the point groups not included are H I and RI I; they need a separate complete analysis, and due to this, they are postponed for a future work. Table 1. Estimated order of magnitude of the high-order tensors for the nonlinear effects of AlN and PZT − 5H (4 mm and 6 mm point groups respectively), using the stress-charge formulation presented.

Nonlinear Effects of Piezoelectric Materials
The nonlinear phenomena of the piezoelectric effect take importance when the material is subject to relatively high electric fields and strong deformations, and its consequences can be classified into two categories. First is the change of the mechanical and electrical properties such as the change of electrical permittivity, elasticity constants, and piezoelectric coupling coefficients. Second is the behavior variation of the physicalsystem response due to the modified material parameters, in particular the arising of the hysteresis behavior, changes in the electromechanical coupling factor, a shift of the resonance frequency, and the modification of the capacitance of the devices.

Variation of Mechanical and Electrical Properties
The change of the electrical permittivity in a piezoelectric material is produced by strong deformations or high temperatures [42] and can be induced by exciting the material with a relatively high electric field, the converse piezoelectric effect producing the deformations needed. This physical behavior can be observed from the nonlinear state equations, since the tensor q jkλ in the D i equation modifies the total polarization, and this can be integrated into a unique term with the electrical permittivity as follows: where e f f ij is the effective electrical permittivity. The change of the elasticity constants is due to exposing the material to relatively high electric fields, which provokes a change in the interatomic electronic forces due to deformations, consequently causing a variation of the stiffness of the material. Furthermore, this phenomenon is included in the state equations through the modification of the total stress induced by the contribution of the deformations and can be formulated as effective elasticity constants as C e f f λµ = C λµ + g λµk E k (25) Finally, due to the power balance of the nonlinear state equations, the variation of the piezoelectric coefficients is a consequence of the imbalance produced by the two last phenomena discussed, where the variation in the transduced power produced by the first nonlinear effect is compensated by the second, then the effective piezoelectric coupling coefficients are where e e f f jλ and e T−e f f iµ must be used in the T λ and D i state equations, respectively.

Change Response of the PhysicalSystem
The literature shows how some piezoelectric devices that are subject to relatively high electric fields have a shift of their resonance frequency; this is produced by the variation of the electrical permittivity and the elasticity constants phenomenon explained before [20,[43][44][45]. Generally, the resonance frequency of a piezoelectric resonator depends on its geometric length and the specific material, often calculated as where λ d is the wavelength of the device, ρ the density of the material, and s the effective elasticity constant (in some cases, it can be called the effective Young's modulus), requiring all subscripts to match the main oscillation mode of the studied device. If we analyze the relative change of the resonance frequency (Equation (27)), we can obtain From Equation (29), it can be noticed how the shift of the resonance frequency is a consequence of the changes of the effective elasticity constants, density, and wavelength of the device, where the last two terms are well known, so they can be neglected [31], because the piezoelectric materials are non-centrosymmetric crystals and the transverse/longitudinal dilatation does not provoke the measured order of magnitude for nonlinear effects.
The variation in the capacitance of the devices is explained through the change in the electrical permittivity phenomenon. Normally, the value of the capacitance of devices that have a dielectric as a piezoelectric material is where C e f f is the effective capacitance of the device, A is the electrodes' contact area, t is the thickness of the piezoelectric, and e f f ii is the effective electrical permittivity over the i axis. Performing an analogous relative variation analysis, then where the last two terms can be neglected, inclusive of the nonlinear effects regime. This is due to the absolute displacement of particles because the order is 0.1Å (theoretical prediction) for the piezoelectric materials under these conditions; hence, dA and dt are not comparable with d e f f ii , since its variation is of the order of thousandths [20]. The calculus for the electromechanical coupling factor k 2 e f f is defined for resonant applications of piezoelectric materials and depends on the oscillation mode of the device, material properties, and specific device geometry. k 2 e f f is a measure of the exchange of power transduced between the mechanical and electrical fields, and for the most common devices, it has an expression of the form [38,46] k 2 for a device with Z-shear oscillation mode and a wave in the X-propagation axis. In Equation (32), the most significant variation, following the discussion above, comes from the effective elasticity constants [31]; hence, when a shift of the resonance frequency occurs, the electromechanical coupling factor increases its value, while the effective elasticity constants decrease. Therefore, for applications where power transduction is the main goal (e.g., energy harvester, microphones, etc.), it can be deduced using Equations (28) to (32) that a negative external electric field increases the performance of the device [13]. The appearance of the hysteresis behavior in the piezoelectric materials as a soft ferroelectric effect is a well-known phenomenon; this is produced by two main causes, the alignment of the dipoles in the unit cells of the material with respect to an external electric field and the change of the domain walls [47][48][49][50]. The change in the domain walls produces a spontaneous strain, inducing additional stress and polarization, and the alignment of the unit cells corresponds to the spontaneous polarization field having a contribution to the strain field. The correspondence between a cyclic electric field and the response of the polarization field and deformation field results in a hysteresis loop and butterfly loop, respectively [51]. The thermodynamic formulation presented only considers the spontaneous strain produced by high electric fields induced due to the inverse piezoelectric effect (last term of the T λ state equations), but it is only one of the theoretical treatments needed for a complete description of the hysteresis behavior.
In the context of all the experimental evidence exposed and discussed, the several nonlinear effects in piezoelectric materials take importance in different regimes. We describe the limit of the formulation presented as a function of the importance of the high-order tensors for their respective regime of operation, where the nonlinear electric contributions take precedence over the mechanical ones [52,53]. Taking as independent physical magnitudes the electrical and deformation field, if the material is subject to excitements of an order of magnitude under 10 6 V/m and 10 −6 , respectively, the linear formulation would be enough. From there, the g λνi and q jkλ tensors must be taken into account, where the r ijk domain makes electric contributions with electric fields above 10 9 V/m, and t λµν is only required starting from the plastic regime. Finally, the hysteresis behavior appears as a soft ferroelectric effect for some specific piezoelectric crystals with excitements of the order above 10 6 V/m and 10 −4 for the electric and deformation fields, respectively. It is necessary to bear in mind that, currently it is not clear what the starting point for the hysteresis behavior for any piezoelectric material is, since this effect belongs to the point group of the material or is induced by very high electric or deformation fields. The last discussion only applies to crystalline piezoelectric materials that are subject to nonlinear perturbative excitements.

Experimental Validation: Simulation
To validate the theoretical development performed in this work, we chose a reference that showed the nonlinear behavior of the piezoelectric devices under a relatively high electric field, since this is the simplest method to induce the nonlinear phenomena. The reference to reproduce is [28], where a solidly mounted resonator (SMR) was fabricated and characterized using AlN as a piezoelectric material; the fabrication details can be found in the reference. Measurements were performed with an Advantest R3767 S-Parameter analyzer, with the DC offset generated by a Keithley K327 and connected through a bias-T, and finally, the data acquisition was performed with the Picoprobe ECP18 GS-200 PP. Therefore, to implement the nonlinear state equations deduced, we can start neglecting the contributions of the t λµν and r ijk tensors, since the operation regime and nonlinear behavior of the device are dominated by the linear description and the tensors g λνi and q jkλ [52,53].
To include the nonlinear state equations within the simulations in an easy way, we chose the following formulation: where we used Equations (25) and (26), since this implementation included the power balance between the physics magnitudes of interest (T λ and D i ) within the FEM simulator. Generally, the FEM simulators allow us to set the electrical relative permittivity as an input parameter, so we rewrite where 0 is the vacuum electrical permittivity, r ij is the relative permittivity in the linear regime, r−e f f ij is the effective relative permittivity, and the last term is defined as The SMR devices have the main oscillation mode, which confines the mechanical waves within the device; based on this, the algebraic tensor development of Equation (33) results in the only components of g λνi and q jkλ that must be taken into account to be g 333 , q r 331 , and q r 333 . The symmetry structure taken from Tables 3 and 4 was the 6 mm one, since the piezoelectric material was AlN. The values obtained for the high-order tensors from the simulations were g 333 = −80N/Vm, and q r 331 = q r 333 = −120 (36) The SMR device was powered by an S-Parameter analyzer with a DC bias added with a bias-T through the signal probe of the RF microprobes. To reproduce the experimental setup, the simulated device was connected to an RF source with a DC voltage overlap, to calculate the whole interest frequency spectrum as a function of the DC bias. Figure 1a shows a transversal cut of the device simulated. In Figure 1b, the impedance of the device simulated for different DC biases is presented, and there, we can observe how the frequency response depends on the external DC electric field (EDEF), since it augments the stiffness of the material when positive voltages are applied, increasing the elasticity constants' values, and consequently, the resonance frequency increases as well; this behavior's prediction was performed by Equation (27). The effective elasticity constant C D−e f f 33 obtained from the measurements and simulations is shown in Figure 2a, where the maximum percent error obtained was 0.79%. The stiffening of the material was proportional to the EDEF due to the negative sign of g 333 ; consequently, the resonance frequency had the same dependency. This can be observed in Figure 2b, where the resonance frequencies measured and simulated, for several values of the EDEF, are presented; there, the maximum percent error obtained was 0.3%. The behavior obtained from the measurements and simulations for the relative effective permittivity is exposed in Figure 3a, where the linear inverse dependence between the EDEF and the permittivity can be observed, as predicted by Equation (24); the maximum percent error obtained was 2.9%. As can be expected, the slopes in Figures 2a and 3a correspond to the values of g 333 and q 331 , having the correct prediction for the trend behavior observed experimentally. Finally, the behavior of the electromechanical coupling factor is shown in Figure 3b, where the predictions of Equations (28) to (32) are corroborated, since the maximum value for k 2 e f f was obtained under negative voltages for the EDEF; this behavior was not reported by the experimental reference, but it was obtained from the simulations. In Table 5 is shown the average and maximum percent errors obtained from the simulations with respect to the measurements; there, the maximum percent error for the effective elasticity constants, effective relative permittivity, and resonance frequencies were 0.79%, 2.9%, and 0.3% respectively. These errors were caused by the difference between the physical material parameters and those used in the simulations; furthermore, the inaccuracy in the extremes of the values of the EDEF was due to the divergence problems that are present in the direct solver of the FEM software; this can be observed mainly in Figures 1b and 3a. Nevertheless, in the scope of the simulations performed, the maximum percent error obtained for any material parameter or physicalsystem parameter was 1.1%; this shows the accuracy of the state equations presented to predict the main nonlinear phenomena of piezoelectric materials through a unified set of state equations, which can be included in FEM simulators easily.   r−e f f 33 obtained from the simulations and measurements of the device fabricated in [28] for an external DC electric field in the range of −2 MV/m to 2 MV/m. (b) Electromechanical coupling factor for the device simulated; this parameter was not reported by [28].

Conclusions
In this work, we presented the nonlinear state equations for piezoelectric materials obtained from first-principles, conserving the power balance exchange between the dependent physical magnitudes T λ and D i and having a unified set of equations that predicts the behavior of the nonlinear phenomena. Furthermore, we showed how we obtained the transformation laws and the symmetry structures for the r ijk , g λµk , and q ijλ tensors, while the calculation procedure was demonstrated with known tensor structures (C λµ , e λk , and ij ). The physical connection and explanation for the nonlinear phenomena experimentally observed in the piezoelectric material were exposed, remarking on the excitement conditions that made each phenomenon appear, where, under an external DC electric field less of than 10 9 V/m, the nonlinear phenomena were dominated by the change in the relative effective permittivity and effective elasticity constants through the g λµk and q ijλ tensors. The elastoelectric effect does not appear in non-piezoelectric materials (g λµk is null), but the electrostrictive effect and nonlinear piezoelectric behavior remained within the material since qijλ was not zero, except for point group 432, where qijλ = 0 and g λµk was not null. A fast methodology for the implementation of the nonlinear state equations in the main FEM simulation software was exposed and demonstrated; this was carried out through the reproduction of an experimental reference, where the main nonlinearities of the piezoelectric effect were measured. The maximum percent errors obtained from the simulations were 0.79%, 2.9%, and 0.3% for the effective elasticity constants, relative effective permittivity, and resonance frequencies. This proved the effectiveness of the nonlinear stress-charge formulation presented, taking into account that the symmetry structure of each high-order tensor was shown (Tables 3, 4   Acknowledgments: A.F. Jaramillo Alvarado thanks Instituto Nacional de Astrofísica, Óptica y Electrónica (INAOE) for its support through the scholarship Beca Colaboración.

Conflicts of Interest:
The authors declare no conflict of interest.